home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / bin / geomutil / math2oogl / BezierPlot.m next >
Text File  |  1993-10-09  |  3KB  |  63 lines
  1. (* Author: Silvio Levy *)
  2.  
  3. BeginPackage["BezierPlot`"]
  4.  
  5. BezierPatch::usage = "\n
  6.   BezierPatch[ 4 by 4 array ] describes a bicubic Bezier patch\n
  7.   to be plotted using OOGL.m, for example.  The elements of the\n
  8.   array are usually 3D vectors, but the function doesn't care."
  9.  
  10. BezierPlot::usage = "\n
  11.   BezierPlot[f,{t,tmin,tmax},{u,umin,umax}] \"plots\" the function f\n
  12.   of two variables t and u, using Bezier patches.  (Usually f returns\n
  13.   a 3D vector.)  BezierPlot returns an array of BezierPatch'es.\n
  14.   Options:\n
  15.     Epsilon is the increment used to compute derivatives (default: 10^-4)\n
  16.     PlotPoints is a number or pair of numbers indicating the fineness\n
  17.       of the subdivision (default: 5)";
  18.  
  19. BezierPlot::baditer = "bad iterator";
  20. BezierPlot::badpp = "requested number of plot points is < 1";
  21.  
  22. Options[BezierPlot] = {Epsilon -> .0001, PlotPoints->5}
  23.   
  24. Begin["`private`"]
  25.  
  26. BezierPlot[f_,{t_,tmin_,tmax_},{u_,umin_,umax_},options___]:=
  27. Block[{it,iu,nt,nu,dt,du,loc,dfdt,dfdt0,dfdu,dfdu0,dboth,all,
  28.       plotPoints=PlotPoints/.{options}/.Options[BezierPlot],
  29.       eps=Epsilon/.{options}/.Options[BezierPlot]},
  30.   If[Length[plotPoints]==2,{nt,nu}=plotPoints,nt=plotPoints;nu=plotPoints];
  31.   If[nx<1 || ny<1,Message[BezierPlot::badpp];Return[{}]];
  32.   dt=N[(tmax-tmin)/nt]; du=N[(umax-umin)/nu];
  33.   If[!(NumberQ[dt]&&dt>0&&NumberQ[du]&&du>0),
  34.     Message[BezierPlot::baditer];Return[{}]];
  35.   loc=Table[f/.{t->tmin+dt it,u->umin+du iu},{iu,0,nu},{it,0,nt}];
  36.   dfdt0=-loc+Table[f/.{t->tmin+dt(it+If[it<nt,1,-1]eps),u->umin+du iu},
  37.                    {iu,0,nu},{it,0,nt}];
  38.   dfdt=MapAt[-1#&,#,-1]& /@ dfdt0/(3 eps);
  39.   dfdu0=-loc+Table[f/.{t->tmin+dt it,u->umin+du(iu+If[iu<nu,1,-1]eps)},
  40.                    {iu,0,nu},{it,0,nt}];
  41.   dfdu=MapAt[-1#&,#,-1]& @ dfdu0/(3 eps);
  42.   dboth=MapAt[-1#&,#,-1]& /@ MapAt[-1#&,#,-1]& @ 
  43.         (-loc-dfdt0-dfdu0+Table[f/.{t->tmin+dt(it+If[it<nt,1,-1]eps),
  44.                                      u->umin+du(iu+If[iu<nu,1,-1]eps)}, 
  45.                                 {iu,0,nu},{it,0,nt}])/(9 eps^2);
  46.   all=Transpose[{loc,dfdt,dfdu,dboth},{3,1,2}]//N;
  47.   Map[
  48.   BezierPatch[
  49.   {{#[[1,1]],#[[1,1]]+#[[1,2]],#[[2,1]]-#[[2,2]],#[[2,1]]},
  50.   {#[[1,1]]+#[[1,3]],#[[1,1]]+#[[1,2]]+#[[1,3]]+#[[1,4]],
  51.     #[[2,1]]-#[[2,2]]+#[[2,3]]-#[[2,4]],#[[2,1]]+#[[2,3]]},
  52.   {#[[3,1]]-#[[3,3]],#[[3,1]]+#[[3,2]]-#[[3,3]]-#[[3,4]],
  53.    #[[4,1]]-#[[4,3]]-#[[4,2]]+#[[4,4]], #[[4,1]]-#[[4,3]]},
  54.   {#[[3,1]],#[[3,1]]+#[[3,2]],#[[4,1]]-#[[4,2]],#[[4,1]]}}]&,
  55.     Transpose[{Drop[#,-1]&/@ Drop[#,-1]& @ all, Rest/@Drop[#,-1]& @ all,
  56.                     Drop[#,-1]&/@ Rest @ all, Rest/@Rest @ all}, {3,1,2}],{2}]
  57. ]
  58.  
  59. Format[x_BezierPatch]:="-BezierPatch-"
  60.  
  61. End[]
  62. EndPackage[]
  63.